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Abstract 

Based  on  time-domain  Boussinesq  equations,  a  comprehensive  model  is  con¬ 
structed  to  simulate  wave  shoaling,  wave  breaking  and  wave  runup  in  coastal  re¬ 
gions.  To  illustrate  the  accuracy  of  various  forms  of  Boussinesq  models  for  large 
effect  of  nonlinearity,  solitary  wave  solutions  corresponding  to  several  height  to 
depth  ratios  are  obtained  numerically.  Comparisons  of  permanent  solutions  are 
made  between  Boussinesq  models  and  other  existing  closed-form  expressions. 
Additional  energy  dissipation  terms  are  included  in  the  Boussinesq  equations  to 
account  for  the  effects  of  wave  breaking  and  bottom  friction,  which  play  an  im¬ 
portant  role  for  wave  transformation  in  surf  zone  and  swash  zone  areas.  Model 
results  model  are  compared  to  data  from  a  laboratory  study  of  random  wave 
breaking  over  a  constant  slope  (Mase  and  Kirby,  1992).  Good  agreement  be¬ 
tween  numerical  results  and  experimental  data  is  found. 

Introduction 

Accurate  prediction  of  wave  transformation  from  deep  to  shallow  water  is  im¬ 
portant  to  the  understanding  of  coastal  processes.  As  waves  propagate  from  the 
deep  ocean  towards  coastal  regions,  the  effect  of  bottom  topography  causes  wave 
height  and  wave  shape  to  change  accordingly.  A  combination  of  wave  refraction- 
diffraction,  reflection  and  nonlinear  interaction  takes  place.  The  increase  of  wave 
height  and  the  steepening  of  wave  crest  eventually  lead  to  wave  breaking,  result¬ 
ing  in  the  generations  of  large  scale  turbulent  motion  and  nearshore  circulation, 
as  well  as  the  movement  of  sediment  along  beaches  and  coastlines. 

With  the  inclusions  of  lowest  order  effects  of  nonlinearity  and  frequency 
dispersion,  Boussinesq  equations  provide  a  sound  and  increasingly  well-tested 
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basis  for  the  simulation  of  wave  propagation  in  coastal  regions.  By  using  depth- 
averaged  velocity  as  a  dependent  variable,  Peregrine  (1967)  derived  the  Boussi- 
nesq  equations  for  variable  water  depth.  Numerical  models  based  on  Peregrine’s 
equations  or  equivalent  formulations  have  been  shown  to  give  predictions  which 
compare  quite  well  with  field  data  (Elgar  and  Guza,  1985)  and  laboratory  data 
(Goring,  1978;  Liu,  Yoon  and  Kirby,  1985;  Rygg,  1988),  when  applied  within 
their  range  of  validity. 

In  recent  years,  efforts  have  been  made  to  extend  the  validity  range  of  stan¬ 
dard  Boussinesq  equations  to  intermediate-depth  areas.  Madsen  et  al.  (1991) 
improved  the  linearized  model  by  introducing  expressions  in  the  equations  which 
are  formally  equivalent  to  zero  within  the  accuracy  of  the  model,  thus  obtain¬ 
ing  a  rearrangement  of  higher-order  terms  in  the  momentum  equations.  Nwogu 
(1993)  used  the  velocity  at  a  certain  depth  as  a  dependent  variable  and  pursued 
a  consistent  derivation  of  the  governing  equations  using  this  non-standard  de¬ 
pendent  variable.  The  resulting  dispersion  relations  of  Madsen  et  al.  (1991)  and 
Nwogu(1993)  are  formally  equivalent  and  are  much  closer  to  the  exact  solution 
in  intermediate  water  depths  than  are  the  standard  Boussinesq  equations. 

By  using  the  velocity  at  a  certain  water  depth  as  a  dependent  variable 
(Nwogu,  1993)  and  making  no  assumption  for  small  nonlinearity,  Wei  et  al. 
(1995)  derived  fully  nonlinear  Boussinesq  equations  which  further  extend  the 
range  of  validity  for  Boussinesq  models.  In  addition  to  obtaining  the  same  dis¬ 
persion  relation  in  intermediate  water  depths  as  those  of  Madsen  et  al.  (1991) 
and  Nwogu  (1993),  the  fully  nonlinear  Boussinesq  equations  provide  an  improved 
theory  for  second  and  third  order  nonlinear  interactions,  and  can  be  applied  to 
simulate  wave  propagation  prior  to  wave  breaking  where  nonlinearity  is  expected 
to  be  large.  Numerical  examples  given  by  Wei  et  al.  (1995)  have  shown  that 
the  fully  nonlinear  Boussinesq  model  predict  wave  heights,  wave  shapes  and 
wave  celerities  much  more  accurately  compared  to  those  from  weakly  nonlinear 
extended  Boussinesq  equations. 

To  extend  Boussinesq  models  to  surf  zone  and  swash  regions  where  the  ef¬ 
fects  of  wave  breaking  and  bottom  friction  are  important,  additional  energy 
dissipation  terms  are  included  to  the  model  equations.  Observing  the  similarity 
between  broken  wave  propagation  and  hydraulic  jumps,  Heitner  and  Housner 
(1970)  proposed  an  eddy  viscosity  model  to  dissipate  energy  for  breaking  waves. 
Energy  loss  is  limited  to  the  front  face  of  waves  where  the  change  of  wave  prop¬ 
erties  exceeds  a  certain  criteria.  Zelt  (1991)  implemented  the  eddy  viscosity 
formula  in  a  Lagrangian  Boussinesq  model  to  simulate  solitary  wave  breaking 
and  runup.  Good  agreements  between  numerical  results  and  experimental  data 
were  obtained. 

In  this  study,  we  are  constructing  a  coastal  processes  model  based  on  the 
extended  Boussinesq  equations.  With  the  improved  dispersion  relation  in  in¬ 
termediate  water  depth,  the  model  can  be  applied  in  relatively  deeper  water. 
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Energy  dissipation  terms  similar  to  those  used  by  Heitner  and  Housner  (1970) 
and  by  Zelt  (1991)  are  added  to  the  Boussinesq  model  to  simulate  wave  propaga¬ 
tion  in  surf  zone  and  swash  zone  areas.  The  same  high-order  predictor-corrector 
finite  difference  scheme  as  that  in  Wei  et  al.  (1995)  is  used  to  discretize  the 
model  equations  and  to  obtain  corresponding  solutions. 

The  Boussinesq  equations  derived  by  Nwogu  (1993)  and  by  Wei  et  al.  (1995) 
are  similar,  with  extra  higher  order  nonlinear  dispersive  terms  including  in  the 
latter  set  of  equations.  To  illustrate  the  importance  of  these  terms,  we  apply 
the  two  models  to  obtain  solitary  wave  solutions  for  several  cases  of  large  non¬ 
linearity.  Results  with  and  without  high  order  nonlinear  dispersive  terms  are 
compared  with  existing  closed-form  solutions.  To  demonstrate  the  effect  of  en¬ 
ergy  dissipation  terms  included  in  the  model  equations,  we  apply  to  laboratory 
study  of  random  wave  breaking  over  a  slope  (Mase  and  Kirby,  1992).  Compar¬ 
isons  of  numerical  results  and  experimental  data  for  wave  heights,  wave  phase 
and  third  moment  wave  statistics  will  be  provided. 


Model  Description 

The  fully  nonlinear  Boussinesq  equations  of  Wei  et  al.  (1995)  are  derived 
in  a  similar  way  as  that  of  Nwogu  (1993)  by  using  the  velocity  at  a  certain 
water  depth  as  a  dependent  variable.  In  the  derivation,  however,  no  assumption 
was  made  about  the  relative  size  of  nonlinear  effects.  The  resulting  equations 
therefore  include  additional  high  order  nonlinear  dispersive  terms.  The  fully 
nonlinear  Boussinesq  equations  in  dimensionless  form  are 


Vt  +  V  •  |  (h  +  6t}) 


ua  + 
+ 


f  -  |( fi- 2  -  h6rj  +  (^)2))  V(V  •  ua) 

f  [za  +  \{h  -  6 q))  V(V  •  W] }  =  0  (1) 

uat  +  £(ua  •  V)ua  +  V?/  +  y2  '  u°t )  +  zaV[V  •  (feual)]J 

+  5//2V  {  (za  -  SV)(ua  •  V)  [V  •  (/m„)]  -  77  [^r/V  •  uat  +  V  •  (/mat)] 

+  \  \zl-  M2]  K  •  V)(V  •  ua)  +  |  [V  •  (feu*)  +  6r,V  ■  ua]2}  =  0  (2) 


where  7/  is  the  surface  elevation,  h  the  water  depth,  u0  the  horizontal  veloc¬ 
ity  vector  at  the  water  depth  of  z  =  zah ,  V  the  horizontal  gradient  operator, 
subscript  t  the  partial  derivative  in  time,  6  the  ratio  of  typical  wave  height  to 
typical  water  depth,  //2  the  product  square  of  typical  wavenumber  and  typical 
water  depth.  The  two  dimensionless  parameters  S  and  y?  represent  the  effects  of 
nonlinearity  and  frequency  dispersion,  respectively.  For  the  case  of  small  non¬ 
linearity  where  terms  with  order  of  6y2  and  higher  may  be  neglected,  the  above 
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equations  reduce  to  Nwogu’s  extended  Boussinesq  equations.  As  will  be  shown 
in  the  next  section,  however,  it  is  these  high  order  nonlinear  dispersive  terms 
that  provide  accurate  predictions  for  wave  heights  and  wave  shapes  of  solitary 
wave  solutions. 

In  addition  to  the  necessary  corrections  for  waves  with  strong  nonlinearity, 
the  fully  nonlinear  Boussinesq  equations  also  provide  a  physically  correct  condi¬ 
tion  for  mass  balance  at  the  shoreline.  As  shown  in  the  continuity  equation  (1), 
all  the  mass  flux  terms  have  a  common  factor  /i-f  5y.  It  is  thus  clear  that  the  mass 
flux  at  the  shoreline  will  go  to  zero  when  the  total  water  depth  h  +  8r]  becomes 
zero.  This  result  is  expected  on  physical  grounds  and  appears  in  the  nonlinear 
shallow  water  equations  and  in  Boussinesq  equations  where  the  depth- averaged 
velocity  is  the  dependent  variable.  However,  this  condition  is  not  automatically 
satisfied  by  Nwogu’s  or  other  weakly  nonlinear  Boussinesq  model  based  on  a 
velocity  other  than  the  depth-averaged  value,  making  the  application  of  these 
models  problematic  at  the  shoreline. 

Following  the  approach  by  Ileitner  and  Ilousner  (1970)  and  by  Zelt  (1991), 
the  overall  effect  of  wave  breaking  is  simulated  by  adding  the  following  energy 
dissipation  term  into  the  right  hand  side  of  equation  (2) 

Fbreak  —  [(^x)i  T  (t'&Wj^y,  (t/bVx')x  H*  (i^i >Uy)y]  (3) 


where  u  and  v  are  the  x  and  y  component  of  ua.  The  notation  Vb  is  the  eddy 
viscosity  defined  as 

vh  =  -Ba2h2S7  •  ua  (4) 


where  B  is  a  coefficient  related  to  the  local  property  of  the  waves  and  the 
corresponding  critical  value  for  wave  breaking  to  take  place.  The  coefficient 
a  is  the  mixing  length  parameter,  whose  value  is  determined  empirically.  In  the 
computation  which  will  be  shown  below,  we  use  it  as  a  constant,  i.e.  a  =  2. 
Zelt  (1991)  defined  the  critical  value  of  velocity  gradient  as  it*  =  -0.3 yjgjh, 
and  coefficient  B  is  given  by 


if  V  •  u  <  2u*x 
if  2u*  <  V  •  u  <  u* 
if  V  •  u  >  u* 


(5) 


In  swash  zone  areas  where  water  moves  up  and  down  along  the  beach  face, 
the  effect  of  bottom  friction  is  no  longer  small.  However,  exact  expressions  for 
friction  force  are  difficult  to  obtain  and  parameterization  formula  has  to  be  used 
instead.  Here,  we  use  Chevy’s  formula  to  compute  bottom  friction,  i.e.,  adding 
the  following  expression  to  the  right  hand  side  of  equation  (2) 


i?  -  U«KI 

r  bottom-  C2^h  +  S^ 


(6) 
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where  C/  is  the  Chevy’s  dimensional  coefficient  which  is  a  function  of  bottom 
roughness  and  the  magnitude  of  velocity  u0.  For  simplicity,  a  constant  coefficient 
is  used  in  the  model.  To  overcome  the  difficulty  of  defining  physical  variables 
in  the  physically  dry  grids  in  a  Eulerian  system,  a  minimum  thickness  of  water 
is  maintained  for  those  dry  grids.  This  thin  layer  of  water  remains  almost 
motionless  due  to  the  balance  between  gravitational  force  and  bottom  friction. 

A  high-order  numerical  scheme  is  utilized  to  obtain  solutions  to  the  model 
equations  described  above.  The  scheme  is  similar  to  that  of  Wei  et  al.  (1995), 
wjth  the  inclusion  of  extra  energy  dissipation  terms  in  the  governing  equations. 
Standard  five-point  and  three-point  finite  difference  schemes  are  used  to  dis¬ 
cretize  the  first-order  spatial  derivative  terms  and  other  higher  order  terms,  re¬ 
spectively.  For  time-stepping,  we  use  the  fourth-order  Adams-Bashforth-  Moul¬ 
ton  predictor-corrector  scheme.  Detailed  formulations  can  be  found  in  Wei  et 
al.  (1995). 

Results  and  Comparisons 

To  demonstrate  the  importance  of  high  order  nonlinear  dispersion  terms  in 
the  fully  nonlinear  Boussinesq  equations,  we  apply  Boussinesq  models  with  and 
without  those  terms  to  obtain  solitary  wave  solutions  for  large  effects  of  nonlin¬ 
earity.  These  permanent  wave  solutions  are  then  compared  with  other  existing 
solutions  given  by  Tanaka  (1986)  and  by  Seabra-Santos  et  al.  (1987).  The  ap¬ 
proximate  analytical  solution  of  solitary  waves  derived  by  Wei  and  Kirby  (1995) 
is  not  valid  for  this  case  due  to  the  assumption  of  small  nonlinearity  used  in  the 
derivation.  As  far  as  we  know,  there  is  no  closed  form  solutions  for  solitary  wave 
for  the  extended  Boussinesq  equations.  Therefore,  numerical  experiments  axe 
performed  here  to  obtain  solitary  waves  corresponding  to  Boussinesq  equations. 

Using  the  approximate  expressions  of  Wei  and  Kirby  (1995)  as  inputs,  we 
ran  both  models  over  a  long  distance  with  constant  depth.  At  the  beginning  of 
these  computations,  wave  heights  and  wave  shapes  changed  constantly,  and  small 
oscillatory  tails  developed  behind  the  main  waves  as  they  propagated  forward. 
After  running  the  models  for  a  long  time,  however,  the  changes  of  form  became 
negligible  and  wave  shapes  stabilized,  indicating  that  a  numerical  permanent- 
form  solitary  wave  solution  corresponding  to  each  of  the  Boussinesq  models  was 
obtained.  Due  to  the  discrepancy  between  initial  and  stabilized  wave  forms, 
several  runs  of  the  models  were  required  to  obtain  a  solitary  wave  with  desired 
height. 

Closed-form  solutions  of  permanent  waves  corresponding  to  different  wave 
equations  have  been  investigated  in  the  past.  Problems  for  wave  propagation 
in  inviscid  and  impressible  fluid  are  governed  by  the  Laplace  equation  inside 
the  fluid  and  associate  surface  and  bottom  boundary  conditions.  By  assuming 
small  effect  of  dispersion  and/or  nonlinearity,  approximate  forms  of  equations  for 
wave  propagation  are  derived  to  reduce  the  dependence  of  vertical  space.  Shallow 
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water  equations,  Boussinesq  equations  and  the  equations  derived  by  Serre  (1953) 
all  belong  to  approximate  equations.  Permanent  solutions  corresponding  to  the 
original  wave  problem  were  obtained  by  Tanaka  (1986),  who  used  conformal 
mapping  to  transform  the  curved  surface  elevation  into  a  straight  line.  Though 
numerical  iteration  is  required  to  solve  the  closed-  form  expression  for  given  wave 
height,  Tanaka’s  solutions  are  the  most  accurate  to  compare  with.  For  long  waves 
propagating  over  constant  depth,  Serre  (1953)  derived  a  set  of  equations  which 
include  more  high  order  nonlinear  dispersion  terms  than  standard  Boussinesq 
equations,  but  not  as  many  as  the  fully  nonlinear  Boussinesq  equations  of  Wei 
et  al.  (1995).  As  shown  in  the  paper  by  Seabra-Santos  et  al.  (1987),  Serre’s 
equations  has  closed-form  solitary  wave  solutions. 

Figure  1  shows  the  comparisons  between  numerical  solutions  of  solitary  waves 
by  extended  Boussinesq  models  and  the  closed-form  solutions  of  Tanaka  (1986), 
and  the  solutions  of  Serre’s  equations.  The  results  from  the  fully  nonlinear 
Boussinesq  model  match  very  well  Tanaka’s  solution  for  all  three  height  to  depth 
ratios.  The  permanent  solutions  corresponding  to  the  extended  Boussinesq  equa¬ 
tions  without  high  order  nonlinear  dispersive  terms  and  to  Serre’s  equations,  on 
the  other  hand,  predicted  either  narrower  or  wider  solitary  wave  shapes.  The 
results  imply  that  all  the  high  order  nonlinear  dispersive  terms  are  important 
for  simulating  the  propagation  of  waves  with  strong  nonlinearity.  Wei  et  al. 
(1995)  applied  the  Boussinesq  models  to  study  solitary  wave  shoaling  over  dif¬ 
ferent  slopes  and  undular  bore  propagation  over  constant  depth.  Results  of 
wave  height,  wave  shape  and  wave  celerity  from  the  fully  nonlinear  Boussinesq 
models  are  much  better  than  those  from  weakly  nonlinear  forms  of  Boussinesq 
equations,  especially  near  wave  breaking  regions  where  the  effect  of  nonlinearity 
is  no  long  small. 

By  including  energy  dissipation  terms  into  the  extended  Boussinesq  equa¬ 
tions,  we  could  further  extend  the  models  to  surf  zone  and  swash  zone  areas, 
where  the  efFects  of  wave  breaking  and  bottom  friction  are  substantial.  Mase 
and  Kirby  (1992)  conducted  a  laboratory  study  of  random  wave  breaking  over 
an  impermeable  beach  with  constant  slope.  Two  sets  of  random  wave  with  peak 
frequencies  0.75 Hz  (test  one)  and  \S)Hz  (test  two)  were  generated  by  a  piston 
wavemaker.  Random  waves  propagated  over  a  constant  depth  of  47.0  cm  and 
then  over  a  constant  slope  of  1  /20.  Time  series  of  surface  elevation  were  collected 
at  12  gage  locations  along  the  slope,  which  will  serve  as  a  verification  for  the 
comprehensive  numerical  model.  The  corresponding  water  depths  for  the  wave 
gauges  are:  47  cm,  35  cm,  30  cm,  25  cm,  20  cm,  17.5  cm,  15.0  cm,  12.5  cm, 
10.0  cm,  7.5  cm,  5.0  cm,  2.5  cm.  Wave  breaking  was  observed  to  start  at  the 
water  depth  of  17.5  cm. 

The  typical  kh  value  (k  is  the  wavenumber  corresponding  to  the  peak  fre¬ 
quency  and  h  the  water  depth)  for  test  two  is  close  to  2,  making  it  impossible 
to  employ  the  standard  Boussinesq  model  due  to  the  dispersion  errors.  The 


x/h0 


•  Figure  1:  Comparison  of  solitary  wave  shapes  for  S  =  0.4, 0.6, 0.8.  Tanaka’s 

solution  ( - );  fully  nonlinear  Boussinesq  model  Boussinesq  model 

without  high-order  nonlinear  dispersive  terms  ( — );  Closed-form  solutions  of 
Serre’s  equations  ( . ) 
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extended  Boussinesq  models,  on  the  other  hand,  could  be  applied  here  to  simu¬ 
late  wave  transformation  due  to  improved  dispersion  relation.  In  the  numerical 
computation,  experimental  data  from  gage  1  were  used  as  input  to  the  model. 
Figure  2  shows  the  comparison  for  surface  elevation  in  the  shallow  regions  where 
wave  breaking  becomes  dominant,  for  a  duration  of  30  seconds.  Except  for  small 
discrepancy  in  wave  phase,  the  model  predicts  the  change  of  wave  form  and  the 
decrease  of  wave  height  due  to  the  effects  of  wave  shoaling  and  breaking.  The 
results  indicate  that  the  model  is  capable  of  simulating  spectral  waves  whose  kh 
values  vary  over  a  large  range. 


time  (sec) 


Figure  2:  Comparison  of  surface  elevation  at  different  water  depths.  Solid  lines 
is  experimental  data  (Mase  and  Kirby,  1992).  Dashed  line  is  numerical  result. 
Wave  breaking  indicated  by  blips  below  the  water  surface  traces. 

To  further  compare  the  wave  properties  between  the  model  and  data,  we 
ran  the  model  for  the  entire  data  set.  Similar  to  those  in  Figure  2,  time  series 
comparisons  are  consistently  good.  Wave  statistics  were  then  computed  from 
the  resulting  time  series.  Figure  3  shows  the  comparisons  of  wave  skewness  and 
asymmetry  for  all  gage  locations.  Except  for  the  last  gage  at  depth  of  2.5  cm 
where  skewness  from  numerical  model  underpredicted  the  data,  good  agreement 
were  found  between  data  and  model.  Since  the  rate  of  sediment  transport  is 


closely  related  to  the  third  moment  of  wave  statistics,  the  model  has  potential 
coastal  engineering  application  for  estimating  sand  movement  along  beaches  and 
coastlines. 


Figure  3:  Comparison  of  skewness  (o)  and  asymmetry  (*)  at  different  water 
depths.  Solid  lines  is  experimental  data  (Mase  and  Kirby,  1992).  Dashed  line  is 
numerical  result. 

Conclusion 

The  fully  nonlinear  Boussinesq  equations  have  been  shown  to  be  capable  of 
predicting  permanent  wave  solutions  more  accurately  than  the  extended  Boussi¬ 
nesq  and  Serre’s  equations,  indicating  that  the  high  order  nonlinear  dispersion 
terms  in  the  equations  are  important  for  simulating  wave  transformation  with 
strong  nonlinear  interaction.  By  including  energy  dissipation  terms  in  the  model 
equations,  the  effects  of  wave  breaking  and  bottom  friction  can  be  simulated. 
Comparison  of  numerical  results  with  laboratory  data  for  random  wave  break¬ 
ing  over  a  slope  are  quite  good.  Agreements  have  been  found  not  only  for  wave 
heights  and  wave  shapes,  but  also  for  skewness  and  asymmetry.  This  results  is 
significant  since  sediment  transport  calculation  depends  on  the  third  moment  of 
wave  statistics.  Future  work  will  be  to  extend  the  model  to  2-D  cases  and  to 
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obtain  wave-induced  currents  by  time  averaging  method. 
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